<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of LatticeSimulation</title>
  <meta name="keywords" content="LatticeSimulation">
  <meta name="description" content="Lattice Simulation Model main routine">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt;  <a href="index.html">Code</a> &gt; LatticeSimulation.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for Code&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->

<h1>LatticeSimulation
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>Lattice Simulation Model main routine</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>function [s,cs,LL,F,Energy]=LatticeSimulation(varargin) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment">Lattice Simulation Model main routine
This function is the main method used to run simulations
This is the main function/routine used to start and compute the simulation. 
Given the temperature,  composition and size of lattice, and the number of
steps, 
this function will pseudo-kinetically step through lattice states
function [s,cs,LL,F,Energy]=Ising2D(kbT,phi,s,Steps)</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
<li><a href="Ising2Doutput.html" class="code" title="">Ising2Doutput</a>	</li><li><a href="LatticeCorrelation.html" class="code" title="function [g]=LatticeCorrelation(s,n,s2)">LatticeCorrelation</a>	Lattice Correlation</li><li><a href="classifypairs2.html" class="code" title="function [pairs,currentEnergy]=classifypairs2(state,energyMatrix,moves,ne,pts)">classifypairs2</a>	Classify the bond pair by calculating change in energy if swapped</li><li><a href="clusterCountEHK2.html" class="code" title="function [c LL F pts]=clusterCountEHK2(s,value)">clusterCountEHK2</a>	Find the clusters of connected'values' in the given lattic state s</li><li><a href="computeInteractionParameters.html" class="code" title="function [ et ] = computeInteractionParameters( chiSA,chiSC,epsilon,gamma )">computeInteractionParameters</a>	computeInteractionParameters calcuate the individual interaction between</li><li><a href="createInitialState.html" class="code" title="function s=createInitialState(phi,size)">createInitialState</a>	createInitialState creates the initial lattice state as a matrix</li><li><a href="neighbors2d.html" class="code" title="function n=neighbors2d(i,j,rows,columns)">neighbors2d</a>	neighbors2d creates a matrix whose entries return a vector of nearest</li><li><a href="plotLattice.html" class="code" title="function a=plotLattice(sites,varargin)">plotLattice</a>	Plot Lattice</li><li><a href="properLabel.html" class="code" title="function r=properLabel(Sn,LL)">properLabel</a>	</li></ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [s,cs,LL,F,Energy]=LatticeSimulation(varargin)</a>
0002 <span class="comment">%Lattice Simulation Model main routine</span>
0003 <span class="comment">%This function is the main method used to run simulations</span>
0004 <span class="comment">%This is the main function/routine used to start and compute the simulation.</span>
0005 <span class="comment">%Given the temperature,  composition and size of lattice, and the number of</span>
0006 <span class="comment">%steps,</span>
0007 <span class="comment">%this function will pseudo-kinetically step through lattice states</span>
0008 <span class="comment">%function [s,cs,LL,F,Energy]=Ising2D(kbT,phi,s,Steps)</span>
0009 
0010 
0011 <span class="comment">%LatticeSimulation(kbT,phi,s,Steps,chiAB,chiAC,epsilon,gamma)</span>
0012 
0013 <span class="comment">% Create an instance of the inputParser class.</span>
0014 p = inputParser;
0015 p.addParamValue(<span class="string">'kbT'</span>,1,@(x)isnumeric(x) &amp;&amp; isscalar(x) &amp;&amp; x&gt;0);
0016 p.addParamValue(<span class="string">'phi'</span>,.5,@(x)isnumeric(x) &amp;&amp; isscalar(x) &amp;&amp; x&gt;=0 &amp;&amp;x&lt;=1);
0017 p.addParamValue(<span class="string">'s'</span>,[64 64]);
0018 p.addParamValue(<span class="string">'Steps'</span>,1000,@(x)isnumeric(x) &amp;&amp; isscalar(x));
0019 p.addParamValue(<span class="string">'chiAB'</span>,-0.1,@(x)isnumeric(x) &amp;&amp; isscalar(x));
0020 p.addParamValue(<span class="string">'chiAC'</span>,0.6,@(x)isnumeric(x) &amp;&amp; isscalar(x));
0021 p.addParamValue(<span class="string">'epsilon'</span>,-1,@(x)isnumeric(x) &amp;&amp; isscalar(x));
0022 p.addParamValue(<span class="string">'gamma'</span>,1.1,@(x)isnumeric(x) &amp;&amp; isscalar(x));
0023 
0024 <span class="comment">%Output Parameters</span>
0025 p.addParamValue(<span class="string">'Sps'</span>,100,@(x)isnumeric(x) &amp;&amp; isscalar(x)); <span class="comment">%steps per save of data</span>
0026 p.addParamValue(<span class="string">'Spp'</span>,3000,@(x)isnumeric(x) &amp;&amp; isscalar(x)); <span class="comment">%steps per plot</span>
0027 p.addParamValue(<span class="string">'Spf'</span>,1500,@(x)isnumeric(x) &amp;&amp; isscalar(x)); <span class="comment">%steps per frame used in creating movies</span>
0028 p.addParamValue(<span class="string">'showplot'</span>,true,@(x)islogical(x) &amp;&amp; isscaler(x));
0029 p.addParamValue(<span class="string">'saveplot'</span>,false,@(x)islogical(x) &amp;&amp; isscaler(x));
0030 p.addParamValue(<span class="string">'saveFrames'</span>,true,@(x)islogical(x) &amp;&amp; isscaler(x));
0031 p.addParamValue(<span class="string">'showMoves'</span>,false,@(x)islogical(x) &amp;&amp; isscaler(x));
0032 
0033 p.parse(varargin{:});
0034 <span class="comment">% Show the value of a specific argument.</span>
0035 params=p.Results;
0036 clear p;
0037 
0038 <span class="comment">%Random Number Generator Setup</span>
0039 <span class="comment">%rand('twister',sum(100*clock))</span>
0040 <span class="comment">%randn('state',sum(100*clock))</span>
0041 <span class="comment">%rand('twister',10)</span>
0042 <span class="comment">%randn('state',5)</span>
0043 
0044     
0045 et=<a href="computeInteractionParameters.html" class="code" title="function [ et ] = computeInteractionParameters( chiSA,chiSC,epsilon,gamma )">computeInteractionParameters</a>(params.chiAB,params.chiAC,params.epsilon,params.gamma);
0046 
0047 NUMBER_NEIGHBORS=8;
0048 
0049 states=[1 2 3 4 5];
0050 <span class="comment">%1--solvent</span>
0051 <span class="comment">%2 --amorphous</span>
0052 <span class="comment">%3 --Crystal</span>
0053 <span class="comment">%4 --Wall</span>
0054 <span class="comment">%5 --Water</span>
0055 
0056 <span class="comment">%Set up Default Arguments</span>
0057 
0058 
0059 disp(<span class="string">'Setting up Initial state'</span>)
0060 
0061 
0062     <span class="comment">%If we input the size of the grid</span>
0063     <span class="keyword">if</span> length(params.s)==2
0064         rows=params.s(1);
0065         cols=params.s(2);
0066         s=<a href="createInitialState.html" class="code" title="function s=createInitialState(phi,size)">createInitialState</a>(params.phi,[rows cols]);
0067     <span class="keyword">else</span>
0068         [rows,cols]=size(params.s);
0069         s=params.s;
0070     <span class="keyword">end</span>
0071 
0072 
0073 
0074 startTime=datestr(now,<span class="string">'yyyymmmdd-HHMM'</span>);
0075 animate=false;
0076 
0077 <span class="keyword">if</span> params.saveFrames
0078     mkdir([<span class="string">'images/'</span> startTime]);
0079     directory = strcat(<span class="string">'images/'</span>,startTime,<span class="string">'/'</span>); 
0080     filename = [directory <span class="string">'Em_T_'</span> num2str(params.kbT) <span class="string">'.txt'</span>];
0081     save(filename,<span class="string">'et'</span>,<span class="string">'-ASCII'</span>)
0082     filename = [directory, num2str(0,<span class="string">'%07d'</span>),<span class="string">'.png'</span>];
0083     imwrite(s,[1 1 1;0 0 0;0 1 0;.5 .5 .5; 0 0 1],filename)
0084 <span class="keyword">end</span>
0085 
0086 orientation=zeros(size(s));
0087 
0088 checkAbsorb=all(s(1,:)==5);
0089 
0090 N=rows*cols;
0091 
0092 <span class="comment">%Initialize Arrays</span>
0093 Energy=zeros(params.Steps+1,1);
0094 frac2=zeros(params.Steps,1);
0095 frac3=zeros(params.Steps,1);
0096 
0097 <span class="comment">%4 neighbors in 2d</span>
0098 <span class="comment">%8 next nearest neighbors in 2d</span>
0099 neighbors=zeros(NUMBER_NEIGHBORS,rows,cols);
0100 
0101 <span class="comment">%how many 1's do we have</span>
0102 m=sum(sum(s==1));
0103 
0104 
0105 <a href="plotLattice.html" class="code" title="function a=plotLattice(sites,varargin)">plotLattice</a>(s,states)
0106 title(<span class="string">'Initial State'</span>)
0107 
0108 <span class="keyword">if</span> params.saveplot
0109     print(<span class="string">'-r600'</span>,<span class="string">'-depsc'</span>,[filestr <span class="string">'InitialPos'</span>]);
0110     close
0111 <span class="keyword">end</span>
0112 
0113 
0114 <span class="comment">%get the possible energy types</span>
0115 <span class="comment">%d_en=testState([-1 1 2]);</span>
0116 
0117 
0118 <span class="comment">%Create Neighbor lists</span>
0119 <span class="keyword">for</span> i=1:rows
0120     <span class="keyword">for</span> j=1:cols
0121         neighbors(:,i,j)=<a href="neighbors2d.html" class="code" title="function n=neighbors2d(i,j,rows,columns)">neighbors2d</a>(i,j,rows,cols);
0122     <span class="keyword">end</span>
0123 <span class="keyword">end</span>
0124 <span class="comment">%Index of Neighbors</span>
0125 <span class="comment">%xxx516xxx</span>
0126 <span class="comment">%xxx3o4xxx</span>
0127 <span class="comment">%xxx728xxx</span>
0128 
0129 <span class="comment">%Not yet fully implemented</span>
0130 <span class="comment">%need to work out details to classify possible moves</span>
0131 swapMoves=[2 4 5 6];
0132 flipMoves=[2];
0133 nMoves=length(swapMoves)+length(flipMoves);
0134 
0135 <span class="comment">%Classify pair types;</span>
0136 <span class="keyword">global</span> pairsMap
0137 pairsMap=-1000*ones(3,3,3,3,3,3,3,3,3,3,3,3);
0138 [px Energy(1)]=<a href="classifypairs2.html" class="code" title="function [pairs,currentEnergy]=classifypairs2(state,energyMatrix,moves,ne,pts)">classifypairs2</a>(s,et,swapMoves,neighbors);
0139 d_en=unique(round(px*1e5)/1e5);
0140 <span class="comment">%Create lists for each scenerio pair</span>
0141  
0142 <span class="keyword">for</span> i=1:length(d_en)
0143     p{i}=(abs(px-d_en(i))&lt;1e-5);
0144     pLength(i)=nnz(p{i});
0145 <span class="keyword">end</span>
0146 
0147 
0148 <span class="comment">%Create energy for each sceneriio pair</span>
0149 <span class="comment">%ignore first value in d_en since it indicates swapping same state</span>
0150 d_energy=[0; d_en(2:end)];
0151 zerogroup=find(d_en==0);
0152 <span class="keyword">if</span> isempty(zerogroup)
0153         zerogroup=length(p)+1;
0154         d_en(zerogroup)=0;
0155         d_energy(zerogroup)=0;
0156         p{zerogroup}=false(size(px));
0157         pLength(zerogroup)=0;
0158         
0159 <span class="keyword">end</span>
0160 d_p=exp(-d_energy./params.kbT);
0161 
0162 cor(:,:,1)=<a href="LatticeCorrelation.html" class="code" title="function [g]=LatticeCorrelation(s,n,s2)">LatticeCorrelation</a>(s,5);
0163 <span class="comment">%s_old=s;</span>
0164 [cs,LL,F]=<a href="clusterCountEHK2.html" class="code" title="function [c LL F pts]=clusterCountEHK2(s,value)">clusterCountEHK2</a>(s,[2 3]);
0165 cs=<a href="properLabel.html" class="code" title="function r=properLabel(Sn,LL)">properLabel</a>(cs,LL);
0166 largest_label=length(LL);
0167 larC(1)=max(LL);
0168 nnc(1)=nnz(LL&gt;5);
0169 fre1(1)=nnz(LL==1);
0170 <span class="comment">%maxLabel=max(labels)+1;</span>
0171 <span class="comment">%clust(1)=sum(csize&gt;10);</span>
0172 <span class="comment">%cave(1)=mean(csize(csize&gt;10));</span>
0173 
0174 disp(<span class="string">'Starting MainLoop'</span>)
0175 <span class="comment">%Number of steps</span>
0176 oldavE=Energy(1);
0177 release=0;
0178 sttt=1;
0179 relt=0;
0180 
0181 <span class="keyword">for</span> t=1:params.Steps
0182     
0183     <span class="comment">%Absorbing top boundry</span>
0184     <span class="keyword">if</span> checkAbsorb
0185     atb=(s(1,:)==2);
0186     <span class="keyword">if</span> any(atb)
0187         sttt=sttt+1;
0188         release(sttt)=release(sttt-1)+nnz(atb);
0189         relt(sttt)=t;
0190         s(1,atb)=5;
0191     <span class="keyword">end</span>
0192     <span class="keyword">end</span>
0193     <span class="comment">%Intermittently save the current state for output</span>
0194     frac2(t)=nnz(s==2);
0195     frac3(t)=nnz(s==3);
0196     <span class="keyword">if</span> mod(t,params.Sps)==0
0197         t
0198         <span class="comment">%LatticeCorrelation(s,0,s_old)</span>
0199         <span class="comment">%avE=mean(Energy(t-Sps+1:t));</span>
0200         <span class="comment">%avE-oldavE</span>
0201         <span class="comment">%oldavE=avE;</span>
0202         [cs,LL,F]=<a href="clusterCountEHK2.html" class="code" title="function [c LL F pts]=clusterCountEHK2(s,value)">clusterCountEHK2</a>(s,[2 3]);
0203         LL23{t/params.Sps+1}=LL;
0204         larC(t/params.Sps+1)=max(LL);
0205         nnc(t/params.Sps+1)=nnz(LL&gt;5);
0206         fre1(t/params.Sps+1)=nnz(LL==1);
0207         [cs,LL,F]=<a href="clusterCountEHK2.html" class="code" title="function [c LL F pts]=clusterCountEHK2(s,value)">clusterCountEHK2</a>(s,[3]);
0208         LL3{t/params.Sps+1}=LL;
0209         <span class="comment">%cor(:,:,t/Sps+1)=LatticeCorrelation(s,5);</span>
0210         <span class="comment">%[c,labels,csize]=clusterCount(s,3);</span>
0211         <span class="comment">%clust(t/Sps+1)=sum(csize&gt;10);</span>
0212         <span class="comment">%cave(t/Sps+1)=mean(csize(csize&gt;10));</span>
0213         
0214     <span class="keyword">end</span>
0215         
0216     <span class="comment">%Running sum of number of pairs * delta energy of switch</span>
0217     <span class="comment">%Skip the first pair category since the swap does not change</span>
0218     q1(1)=pLength(2)*d_p(2);
0219     <span class="keyword">for</span> r=2:length(d_p)-1
0220         q1(r)=q1(r-1)+pLength(r+1)*d_p(r+1);
0221     <span class="keyword">end</span>
0222     
0223     <span class="comment">%Select a random number to determine which category pair to swap</span>
0224        <span class="keyword">if</span> q1(end)==0
0225            <span class="keyword">break</span>;
0226        <span class="keyword">end</span>
0227        temp=rand*q1(end);
0228        st=find(q1&gt;temp,1)+1;
0229        <span class="comment">%Select a random pair from that category</span>
0230        tList=find(p{st});
0231        q=tList(ceil(pLength(st).*rand));
0232        Energy(t+1)=Energy(t)+d_energy(st);
0233        <span class="comment">%Now that q=the pair index lets do the swap and recalculate the</span>
0234        <span class="comment">%list</span>
0235 <span class="comment">%        si=mod(q,N);</span>
0236 <span class="comment">%        if si==0,si=N;end</span>
0237 <span class="comment">%        nd=nmoves(floor((q-1)/N)+1);</span>
0238 <span class="keyword">if</span> q&lt;0
0239     pts=find(<a href="properLabel.html" class="code" title="function r=properLabel(Sn,LL)">properLabel</a>(cs,LL)==-q);
0240     dir=floor(rand*4)+1;
0241     s(pts)=1;
0242     s(neighbors(dir,pts))=2;
0243     newpts=neighbors(dir,pts);
0244     pts=unique([pts;newpts']);
0245     nn=unique(neighbors(:,neighbors(:,pts)));
0246     [cs,LL,F]=<a href="clusterCountEHK2.html" class="code" title="function [c LL F pts]=clusterCountEHK2(s,value)">clusterCountEHK2</a>(s,[2 3]);
0247     cs=<a href="properLabel.html" class="code" title="function r=properLabel(Sn,LL)">properLabel</a>(cs,LL);
0248     largest_label=length(LL);
0249 <span class="keyword">else</span>
0250        [nd,si]=ind2sub(size(px),q);
0251        n_1=neighbors(:,si);
0252        <span class="keyword">if</span> nd&lt;=length(swapMoves)
0253            
0254         nd=swapMoves(nd);
0255         temp=s(si);
0256        
0257         
0258         <span class="keyword">if</span> temp==2 &amp;&amp; s(n_1(nd))==2 <span class="comment">%flip bond</span>
0259             s(si)=3;
0260             s(n_1(nd))=3;
0261             orientation(si)=nd;
0262             orientation(n_1(nd))=nd;
0263         <span class="keyword">else</span>
0264           s(si)=s(n_1(nd));  <span class="comment">%&quot;swap&quot; values</span>
0265         
0266           <span class="comment">%Swap the neighbor of the site that corresponds to the pair</span>
0267           s(n_1(nd))=temp;
0268         <span class="keyword">end</span>
0269 <span class="comment">%        updateClusters</span>
0270         <span class="comment">%get the row and column of the site</span>
0271        <span class="keyword">else</span>
0272            <span class="comment">%we flip states</span>
0273            
0274            nd=1;
0275            <span class="keyword">if</span> s(si)==2
0276                 s(si)=3;
0277            <span class="keyword">else</span>
0278                s(si)=2;
0279                 <span class="keyword">for</span> gi=1:4
0280                     <span class="keyword">if</span> s(neighbors(gi,si))==3
0281                         
0282                         <span class="keyword">if</span> ~any(s(neighbors(1:4,neighbors(gi,si)))==3)
0283                             s(neighbors(gi,si))=2;
0284                         <span class="keyword">end</span>
0285                     <span class="keyword">end</span>
0286                   
0287                 <span class="keyword">end</span>
0288                
0289            <span class="keyword">end</span>
0290           <span class="comment">% n_1=si;</span>
0291            
0292        <span class="keyword">end</span>
0293        
0294        
0295        <span class="comment">%No longer need to look at moves indivdually</span>
0296        
0297        
0298 <span class="comment">%         if showplot &amp;&amp; showMoves</span>
0299 <span class="comment">%             [r c]=ind2sub([rows cols],si);</span>
0300 <span class="comment">%         line('xdata',c,'ydata',-r,'markersize',10,'marker','o', ...</span>
0301 <span class="comment">%     'linestyle','none','markerfacecolor','red','markeredgecolor','black')</span>
0302 <span class="comment">%</span>
0303 <span class="comment">%             [r c]=ind2sub([rows cols],n_1(nd));</span>
0304 <span class="comment">%         line('xdata',c,'ydata',-r,'markersize',10,'marker','o', ...</span>
0305 <span class="comment">%     'linestyle','none','markerfacecolor','red','markeredgecolor','black')</span>
0306 <span class="comment">%</span>
0307 <span class="comment">%         end</span>
0308         
0309         
0310 
0311 
0312         <span class="comment">%for now, reclassify bonds</span>
0313         <span class="comment">%in the future can precalculate changes</span>
0314         
0315         <span class="comment">%find the rows and columns of the neighbors of other pair</span>
0316         
0317         
0318         n_1b=[n_1; neighbors(:,n_1(nd))];
0319         
0320         <span class="comment">%find the neighbors of the neighbors ...</span>
0321         <span class="comment">%these pair energies are</span>
0322         <span class="comment">%afftect by the swap</span>
0323        
0324         nn=unique(neighbors(:,n_1b));
0325         
0326         
0327         <span class="comment">%Loop through the pairs from all the neighbors and</span>
0328         <span class="comment">%reclassify the type</span>
0329         <span class="comment">%update the lists if things change</span>
0330 <span class="keyword">end</span>
0331         [p2,e0]=<a href="classifypairs2.html" class="code" title="function [pairs,currentEnergy]=classifypairs2(state,energyMatrix,moves,ne,pts)">classifypairs2</a>(s,et,swapMoves,neighbors,nn);
0332         
0333 <span class="comment">%         [crap1 EnergyTemp]=classifypairs2(s,et,swapMoves,neighbors);</span>
0334 <span class="comment">%         if abs(EnergyTemp-Energy(t+1))&gt;.01</span>
0335 <span class="comment">%            EnergyTemp</span>
0336 <span class="comment">%            Energy(t+1)</span>
0337 <span class="comment">%         end</span>
0338         <span class="comment">% p{zerogroup}(p{zerogroup}&lt;0)=[];</span>
0339         
0340         <span class="keyword">for</span> kk=1:length(nn)
0341             <span class="keyword">for</span> pp=1:nMoves
0342 
0343                 <span class="keyword">if</span> abs(px(pp,nn(kk))-p2(pp,nn(kk)))&gt;0.001
0344                     oldgroup=find(abs(d_en-px(pp,nn(kk)))&lt;0.001);
0345                     newgroup=find(abs(d_en-p2(pp,nn(kk)))&lt;0.001);
0346                     
0347                     <span class="comment">%Try using logical trackers</span>
0348                     p{oldgroup}(pp,nn(kk))=0;
0349                     pLength(oldgroup)=pLength(oldgroup)-1;
0350                     <span class="comment">%The following search is slow</span>
0351                     
0352 <span class="comment">%                     ind=plistindex(pp,nn(kk));</span>
0353 <span class="comment">%</span>
0354 <span class="comment">%</span>
0355 <span class="comment">%                         plistindex(p{oldgroup}(ind:end))=...</span>
0356 <span class="comment">%                         plistindex(p{oldgroup}(ind:end))-1;</span>
0357 <span class="comment">%</span>
0358 <span class="comment">%                     p{oldgroup}(ind)=[];</span>
0359                     <span class="keyword">if</span> isempty(newgroup)
0360                         newgroup=length(p)+1;
0361                         d_en(newgroup)=p2(pp,nn(kk));
0362                         d_energy(newgroup)=p2(pp,nn(kk));
0363                         d_p=exp(-d_energy./params.kbT);
0364                         <span class="comment">%p{newgroup}=[nMoves*(nn(kk)-1)+pp ];</span>
0365                         p{newgroup}=false(size(px));
0366                         pLength(newgroup)=0;
0367                         
0368                     <span class="keyword">else</span>
0369                         
0370                         <span class="comment">%p{newgroup}=[p{newgroup};nMoves*(nn(kk)-1)+pp ];</span>
0371                         
0372                     <span class="keyword">end</span>
0373                     
0374                         p{newgroup}(pp,nn(kk))=1;
0375                         pLength(newgroup)=pLength(newgroup)+1;
0376                     
0377                     <span class="comment">%plistindex(pp,nn(kk))=length(p{newgroup});</span>
0378                     
0379                 <span class="keyword">end</span>
0380             <span class="keyword">end</span>
0381         <span class="keyword">end</span>
0382         px(:,nn)=p2(:,nn);
0383 <span class="comment">%         for kk=1:length(pI)</span>
0384 <span class="comment">%             if px(pI(kk))~=p2(pI(kk))</span>
0385 <span class="comment">%                 oldgroup=find(d_en==px(pI(kk)));</span>
0386 <span class="comment">%                 newgroup=find(d_en==p2(pI(kk)));</span>
0387 <span class="comment">%                 %The following search is slow</span>
0388 <span class="comment">%                 %should keep track of the index</span>
0389 <span class="comment">%                 p{oldgroup}(p{oldgroup}==pI(kk))=[];</span>
0390 <span class="comment">%                 p{newgroup}=[p{newgroup}; pI(kk)];</span>
0391 <span class="comment">%</span>
0392 <span class="comment">%             end</span>
0393 <span class="comment">%         end</span>
0394 <span class="comment">%         px(pI)=p2(pI);</span>
0395         
0396 <span class="comment">%move clusters</span>
0397         <span class="comment">%[c labs csize cpos]=clusterCount(s,2);</span>
0398         <span class="comment">%bigcI=find(LL&gt;1);</span>
0399    
0400 <span class="comment">%    p{zerogroup}=[p{zerogroup};-bigcI];</span>
0401     
0402         
0403         
0404        <span class="keyword">if</span> params.showplot  &amp;&amp; mod(t,params.Spp)==0
0405        <span class="keyword">if</span> strcmp(get(gcf,<span class="string">'CurrentCharacter'</span>),<span class="string">'p'</span>)    
0406         pause
0407        <span class="keyword">end</span>
0408        
0409        <span class="comment">%subplot(3,3,[1:6])</span>
0410        <a href="plotLattice.html" class="code" title="function a=plotLattice(sites,varargin)">plotLattice</a>(s, states)
0411        <span class="comment">%plotBonds(px,d_en(2:end))</span>
0412        title([<span class="string">'Step'</span> num2str(t,<span class="string">'%05d'</span>)]);
0413        pause(.01)
0414        <span class="comment">%subplot(3,3,7)</span>
0415        <span class="comment">%for i=2:length(p)</span>
0416        <span class="comment">%    ng(i-1)=length(p{i});</span>
0417        <span class="comment">%end</span>
0418        <span class="comment">%bar(d_en(2:end),ng)</span>
0419        <span class="comment">%subplot(3,3,8)</span>
0420        <span class="comment">%plot(0:t,Energy(1:t+1));</span>
0421        <span class="comment">%subplot(3,3,[1:6])</span>
0422        
0423        <span class="keyword">if</span> get(gcf,<span class="string">'CurrentCharacter'</span>)==<span class="string">'p'</span>    
0424         pause
0425        <span class="keyword">end</span>
0426        <span class="keyword">end</span>
0427        <span class="keyword">if</span> params.saveFrames &amp;&amp; mod(t,params.Spf)==0
0428 
0429           <span class="comment">% The next four lines parse and assemble fil</span>
0430 <span class="comment">%</span>
0431          filename = [directory, num2str(t,<span class="string">'%07d'</span>),<span class="string">'.png'</span>];
0432 <span class="comment">%</span>
0433 <span class="comment">%            set(gcf,'PaperPositionMode','auto')</span>
0434             
0435 <span class="comment">%            print(gcf,'-dpng','-r0',filename)</span>
0436             
0437             imwrite(s,[1 1 1;0 0 0;0 1 0;.5 .5 .5;0 0 1],filename)
0438             
0439         <span class="keyword">end</span>
0440 <span class="keyword">end</span>
0441 
0442 <span class="comment">%Plot the outputs</span>
0443 <a href="Ising2Doutput.html" class="code" title="">Ising2Doutput</a>
0444 
0445</pre></div>
<hr><address>Generated on Mon 21-May-2012 09:54:18 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>